Single Cell Investigation

library(ggplot2)
library(ggthemes)
library(reshape2)
library(dplyr)
library(readxl)
library(gridExtra)
theme_set(theme_stata())

sc_eqtl <- read_xlsx("/Users/kathryntsai/OneDrive\ -\ Villanova\ University/College/2018-2019/Summer\ 2019/TFs_eQTLs_Research/RProjects/Project2/IMPACT/NIHMS76345-supplement-Supplementary_table_2_vanderWijst.xlsx", skip=1)

-
/
                                                                                                                                             

ALL SINGLE CELL DATA

gg1 <- melt(sc_eqtl, 
            id.vars=c("SNP", "Gene"), 
            #measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1", "cMono__1", "ncMono__1")
            measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_histogram(alpha=.8, aes(x=-log(value), fill=variable)) + ggtitle("P-values of eQTLs, grouped by cell type - all sc-eQTLs") # possibly take out position="identity" 
gg2 <- melt(sc_eqtl, 
            id.vars=c("SNP", "Gene"), 
            #measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1", "cMono__1", "ncMono__1")
            measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_histogram(position="fill", alpha=.8, aes(x=value, fill=variable)) + ggtitle("Proportions of each cell type at each p-value - all sc-eQTLs") # possibly take out position="identity" 
gg3 <- melt(sc_eqtl, 
            id.vars=c("SNP", "Gene"), 
            measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "cMono__1")
            #measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_bar(alpha=.8, aes(x=value, fill=variable)) + ggtitle("# eQTLs Significant at FDR < 0.05 by Cell Type - all sc-eQTLs") # possibly take out position="identity" 
gg4 <- melt(sc_eqtl, 
            id.vars=c("SNP", "Gene"), 
            measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "cMono__1")
            #measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_bar(position="fill", alpha=.8, aes(x=value, fill=variable)) + ggtitle("Proportions of Cell Types eQTLs Significant at FDR < 0.05") # possibly take out position="identity" 
# No PBMCs
gg5 <- melt(sc_eqtl, 
            id.vars=c("SNP", "Gene"), 
            #measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1", "cMono__1", "ncMono__1")
            measure.vars = c("T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_histogram(alpha=.8, aes(x=-log(value), fill=variable)) + ggtitle("P-values of eQTLs, grouped by cell type - all sc-eQTLs, no PBMCs") # possibly take out position="identity" 
gg6 <- melt(sc_eqtl, 
            id.vars=c("SNP", "Gene"), 
            #measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1", "cMono__1", "ncMono__1")
            measure.vars = c("T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_histogram(position="fill", alpha=.8, aes(x=value, fill=variable)) + ggtitle("Proportions of each cell type at each p-value - all sc-eQTLs, no PBMCs") # possibly take out position="identity" 
gg7 <- melt(sc_eqtl, 
            id.vars=c("SNP", "Gene"), 
            measure.vars = c("T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "cMono__1")
            #measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_bar(alpha=.8, aes(x=value, fill=variable)) + ggtitle("# eQTLs Significant at FDR < 0.05 by Cell Type - all sc-eQTLs, no PBMCs") # possibly take out position="identity" 
gg8 <- melt(sc_eqtl, 
            id.vars=c("SNP", "Gene"), 
            measure.vars = c("T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "cMono__1")
            #measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_bar(position="fill", alpha=.8, aes(x=value, fill=variable)) + ggtitle("Proportions of Cell Types eQTLs Significant at FDR < 0.05 - all sc-eQTLs, no PBMCs") # possibly take out position="identity" 

grid.arrange(gg1, gg2, gg3, gg4, gg5, gg6, gg7, gg8, nrow=2)

How many SNP-gene associations exist in both single cell and bulk datasets?

g <-
  merge(
    sc_eqtl,
    franke_cis_data,
    by = c("SNP", "Gene"),
    all.x = F,
    all.y = F
  )
dim(g)

[1] 317  46

# 317 overlap

Of these overlapping associations, are they more or less cell type specific than the associations that can only be detected in the single cell data?

gg9 <- melt(g, 
           id.vars=c("SNP", "Gene"), 
           #measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1", "cMono__1", "ncMono__1")
           measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
      ) %>% 
      ggplot() + geom_histogram(alpha=.8, aes(x=-log(value), fill=variable)) + ggtitle("P-values of eQTLs, grouped by cell type") # possibly take out position="identity" 
gg10 <- melt(g, 
           id.vars=c("SNP", "Gene"), 
           #measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1", "cMono__1", "ncMono__1")
           measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_histogram(position="fill", alpha=.8, aes(x=value, fill=variable)) + ggtitle("Proportions of each cell type at each p-value") # possibly take out position="identity" 
gg11 <- melt(g, 
            id.vars=c("SNP", "Gene"), 
            measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "cMono__1")
            #measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_bar(alpha=.8, aes(x=value, fill=variable)) + ggtitle("# eQTLs Significant at FDR < 0.05 by Cell Type") # possibly take out position="identity" 
gg12 <- melt(g, 
            id.vars=c("SNP", "Gene"), 
            measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "cMono__1")
            #measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_bar(position="fill", alpha=.8, aes(x=value, fill=variable)) + ggtitle("Proportions of Cell Types eQTLs Significant at FDR < 0.05") # possibly take out position="identity" 
# No PBMCs
gg13 <- melt(g, 
            id.vars=c("SNP", "Gene"), 
            #measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1", "cMono__1", "ncMono__1")
            measure.vars = c("T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_histogram(alpha=.8, aes(x=-log(value), fill=variable)) + ggtitle("P-values of eQTLs, grouped by cell type") # possibly take out position="identity" 
gg14 <- melt(g, 
            id.vars=c("SNP", "Gene"), 
            #measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1", "cMono__1", "ncMono__1")
            measure.vars = c("T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_histogram(position="fill", alpha=.8, aes(x=value, fill=variable)) + ggtitle("Proportions of each cell type at each p-value") # possibly take out position="identity" 
gg15 <- melt(g, 
            id.vars=c("SNP", "Gene"), 
            measure.vars = c("T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "cMono__1")
            #measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_bar(alpha=.8, aes(x=value, fill=variable)) + ggtitle("# eQTLs Significant at FDR < 0.05 by Cell Type") # possibly take out position="identity" 
gg16 <- melt(g, 
            id.vars=c("SNP", "Gene"), 
            measure.vars = c("T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "cMono__1")
            #measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "cMono")
) %>% 
  ggplot() + geom_bar(position="fill", alpha=.8, aes(x=value, fill=variable)) + ggtitle("Proportions of Cell Types eQTLs Significant at FDR < 0.05") # possibly take out position="identity" 

pdf("sc_eqtl_visualization.pdf", height=24, width=50)
grid.arrange(gg9, gg10, gg11, gg12)
grid.arrange(gg13, gg14, gg15, gg16)
dev.off()

Making a Heat Map

https://slowkow.com/notes/pheatmap-tutorial/

library(pheatmap)
library(RColorBrewer)
library(viridis)

gg17 <- melt(sc_eqtl, 
            id.vars=c("SNP", "Gene"), 
            #measure.vars = c("PBMC__1","T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1", "cMono__1", "ncMono__1")
            measure.vars = c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "Mono")
)

mat_col <- data.table(gg17$variable)
mat_colors <- list(group = brewer.pal(length(unique(gg17$variable)), "Set1"))
names(mat_colors$group) <- unique(gg17$variable)

quantile_breaks <- function(xs, n = 10) {
  breaks <- quantile(xs, probs = seq(0, 1, length.out = n))
  breaks[!duplicated(breaks)]
}

mat_breaks <- quantile_breaks(mat, n = 11)


draw_colnames_45 <- function (coln, gaps, ...) {
  coord <- pheatmap:::find_coordinates(length(coln), gaps)
  x     <- coord$coord - 0.5 * coord$size
  res   <- grid::textGrob(
    coln, x = x, y = unit(1, "npc") - unit(3,"bigpts"),
    vjust = 0.75, hjust = 1, rot = 45, gp = grid::gpar(...)
  )
  return(res)
}
assignInNamespace(
  x = "draw_colnames",
  value = "draw_colnames_45",
  ns = asNamespace("pheatmap")
)

gg17

pdf("sc_eqtl_pheatmap.pdf", height=24, width=50)
pheatmap(
  mat               = -log10(sc_eqtl[, c("PBMC","T CD4+", "T CD8+", "NK", "B", "DC", "Mono")]),
  color             = inferno(length(mat_breaks) -1),
  breaks            = mat_breaks,
  border_color      = F,
  cellwidth = 15,
  cellheight= 3,
  legend = TRUE,
  show_colnames     = T,
  show_rownames     = F,
  treeheight_row = 0,
  treeheight_col = 0,
  annotation_names_col    =  data.table(group = unique(gg17$variable)),
  annotation_legend = TRUE,
  #annotation_names_row = data.table(group = paste0(gg17$SNP, "-", gg17$Gene)),
  annotation_colors = mat_colors,
  drop_levels       = TRUE,
  fontsize          = 3,
  main              = "Quantile Color Scale with Log Values"
)
dev.off()

gg19_new <- gg19 %>% 
  # mutate(state=factor(state,levels=rev(sort(unique(state))))) %>%
  mutate(new_value = cut(gg19$value,breaks=8, labels=c("<0.001", 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875))) #%>% 

The `printer` argument is deprecated as of rlang 0.3.0.
This warning is displayed once per session.

gg19_plot + theme_grey(base_size=8) +
  #theme options
  theme(
    #bold font for legend text
    legend.text=element_text(face="bold"),
    #set thickness of axis ticks
    axis.ticks=element_line(size=0.4),
    axis.text.x = element_text(size=1, angle = 45, hjust = 1),
    #remove plot background
    plot.background=element_blank(),
    #remove plot border
    panel.border=element_blank()
  )

png("test_heatmap3.png")#, height = 1000, width = 1280)
y3 <- heatmap.2(-log10(sc_eqtl_plot3), title(main = "cis eQTLs by Cell Type, non-overlap with cis eQTLs", 10), #ColSideColors = map_color_mat_ordered[,2],  labCol = map_color_mat_ordered[,3],
              cex.main = 1, trace = "none", adjCol = 1, adjRow = 0, labCol = cols3,cexRow=1, 
              labRow = rows3, dendrogram = "none", Rowv = F, Colv = F, col = rev(brewer.pal(9,"YlGnBu")),
              cexCol = 1, srtCol = 90, 
              key = T, keysize=2,key.title = "P-value of eQTL for each cell type", key.xlab = "-log10(P-value)",
              density.info = c("none"),margins=c(14, 14),
              sepcolor = "gray", 
              sepwidth=c(0.005,0.005), 
              rowsep=1:ncol(sc_eqtl_plot3),
              colsep=1:nrow(sc_eqtl_plot3))
dev.off()

null device 
          1 

png("test_heatmap.png") # height= 1000, width=1280
y <- heatmap.2(sc_eqtl_plot, title(main = "cis eQTLs by Cell Type", 10), #ColSideColors = map_color_mat_ordered[,2],  labCol = map_color_mat_ordered[,3],
              cex.main = 10, trace = "none", adjCol = 1, adjRow = 0, labCol = cols,cexRow=3, 
              labRow = rows, dendrogram = "none", Rowv = F, Colv = F, #col = rev(brewer.pal(8,"YlGnBu"))
              cexCol = 3, srtCol = 90, 
              key = T, keysize=2,key.title = "P-value of eQTL for each cell type", key.xlab = "P-value",
              density.info = c("none"),margins=c(14, 14),
              sepcolor = "black", 
              sepwidth=c(0.1,0.1), 
              colsep=1:nrow(sc_eqtl_plot),
              rowsep=1:ncol(sc_eqtl_plot)
      )
dev.off()

null device 
          1 

Warning in gzfile(file, "wb") :
  cannot open compressed file '/Users/kathryntsai/OneDrive - Villanova University/College/2018-2019/Summer 2019/TFs_eQTLs_Research/RProjects/Project2/.Rproj.user/shared/notebooks/3CDBF725-single_cell_investigation/1/FEC238FDE19E1822/ce3mtgst2t7bv_t/4fc82b9b2079495bb9e59b2918158cff.snapshot', probable reason 'No such file or directory'
Error in gzfile(file, "wb") : cannot open the connection

g0 <- sc_eqtl[, c("PBMC__1", "T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1")]
cols <- colnames(g0)
sig_sc <- sapply(g0, sub, pattern='\\*', replacement=1)
sig_sc[is.na(sig_sc)] <- 0
g111 <- sig_sc
g111 <- matrix(as.numeric(g111), nrow=nrow(g0), ncol=ncol(g0))
a1 <- data.table(rowSums(g111))
colnames(a1)<- "a"
nrow(a1[which(a1$a == 1)])

[1] 178

g22 <- g[, c("PBMC__1", "T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1")]
cols <- colnames(g22)
g222 <- sapply(g22, sub, pattern='\\*', replacement=1)
g222[is.na(g222)] <- 0
g222 <- matrix(as.numeric(g222), nrow=nrow(g22), ncol=ncol(g22))
a2 <- data.table(rowSums(g222))
colnames(a2)<- "a"
nrow(a2[which(a2$a == 1)])

[1] 137

g33 <- g2[, c("PBMC__1", "T CD4+__1", "T CD8+__1", "NK__1", "B__1", "DC__1", "Mono__1")]
cols <- colnames(g33)
g333 <- sapply(g33, sub, pattern='\\*', replacement=1)
g333[is.na(g333)] <- 0
g333 <- matrix(as.numeric(g333), nrow=nrow(g33), ncol=ncol(g33))
a3 <- data.table(rowSums(g333))
colnames(a3)<- "a"
nrow(a3[which(a3$a == 1)])

[1] 41

Do we even need sc-eQTLs if we can use IMPACT to identify the cell type?

Since Pfizer eQTLs are found in whole blood, can we use IMPACT to infer cell-type-specificity? Since eQTLgen eQTLs are found in peripheral blood, can we use IMPACT to infer cell-type-specificity?

After learning from the annotation analyses above, can we predict if a given SNP will participate in a cis or trans eQTL?

Consider 3 groups of SNPs: only cis, only trans, both cis and trans (with two separate genes, of course)

What are the models of these datasets? Do all include expression PCs? Traditionally, expression PCs don’t take away signal from cis eQTLs but they might from trans eQTLs. Need to know data details, e.g. how many donors? How many samples? How many cells?

proportion of pbmcs are different cells